Drying Induced Hydrophobic Polymer Collapse 
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We have used computer simulation to study the collapse of a hydrophobic chain in water. We find 
that the mechanism of collapse is much like that of a first order phase transition. The evaporation 
of water in the vicinity of the polymer provides the driving force for collapse, and the rate limiting 
step is the nucleation of a sufficiently large vapor bubble. The study is made possible through the 
application of transition path sampling and a coarse grained treatment of liquid water. Relevance 
of our findings to understanding the folding and assembly of proteins is discussed. 
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For nearly a half-century, hydrophobic interactions 
have been considered the primary cause for self assembly 
in soft matter, and a major source of stability in bio- 
physical assembly |^J^. Studying these interactions in 
perhaps their most basic form, we use computer simula- 
tions to demonstrate the mechanism for the collapse of 
a hydrophobic polymer in water. We show that solvent 
fluctuations induce the transition from the extended coil 
to the collapsed globule state, where a vapor bubble of 
sufficient size is formed to envelop and thereby stabilize 
a critical nucleus of hydrophobic units. This mechanism 
is different from that usually considered, where coil to 
globule transitions are attributed to effective interactions 
between pairs of chain segments and a change in sign of 
second virial coefficient |^]. Rather, the mechanism we 
find is evocative of the n-cluster model, where hydropho- 
bic collapse is produced by solvent induced interactions 
between a relatively large cluster of segments Q]. 

As expected from earlier work on the equilibrium the- 
ory of hydrophobicity ||] , we find that the solvent length 
scales pertinent to hydrophobic collapse extend over 
nanometers. We also find that pertinent time scales ex- 
tend beyond nanoseconds. Given these molecularly large 
lengths and times, it is understandable that no work be- 
fore this has provided statistically meaningful computer 
simulations of the process. Our use of a statistical field 
model of water allows us simulate solvent dynamics over 
large length and time scales that would be impractical to 
study with purely atomistic simulation. Spatially com- 
plex small length scale fluctuations are analytically in- 
tegrated out, thus removing the most computationally 
costly features from our simulation. Their integration 
can be performed at the outset because their relaxation is 
relatively fast and their statistics is essentially Gaus- 
sian 0. Only the polymer degrees of freedom and a 
coarse-grained density field remain. The equilibrium the- 
ory for this approach has been detailed in an earlier pa- 
per Here, we use a version of the model that is suit- 
ably generalized for dynamical applications. 

By "small length" we refer to distances smaller than 
I w 0.3 nm. In the absence of any strong perturbation, 
such as those that can occur close to a solute, these small 
length scale fluctuations are the only fluctuations of sig- 



nificance. Larger length scale fluctuations are generally 
insignificant in water at ambient conditions. The liquid 
is relatively cold and incompressible, and spatial correla- 
tions extend over only one or two water molecules. But 
the presence of a hydrophobic polymer can change this 
situation, making large length scale fluctuations impor- 
tant. The liquid lies close to macroscopic vapor-liquid 
equilibrium, and vapor-like behavior is stabilized in the 
vicinity of a sufficiently extended hydrophobic surface , 
a surface formed, for example, by the clustering of hy- 
drophobic groups in a polymer chain. This effect can 
be captured by the behavior of a coarse-grained density 
field , and such a field is conveniently simulated with 
a field of binary numbers jsj, as we discuss now. 



EQUILIBRIUM OF THE SOLVENT MODEL 

Dividing space into a cubic grid of cells, each with side 
length /, this binary field is a dynamical variable given by 
piUi, where pi is a constant given by the bulk liquid den- 
sity, and Hi is dynamical given by = 1 if cell i contains 
"liquid," and = if it contains "vapor" . This coarse- 
grained density is thus the "Ising" field that describes 
phase coexistence and gas-liquid interfaces jl^. The re- 
maining small length scale field, Sp{r), is the difference 
between pirii and the actual density at a point r in cell i. 
It takes on a continuum of values, obeys Gaussian statis- 
tics, and captures the granularity of the solvent [p^ . 

Confining our attention to an ideal hydrophobic solute, 
namely a polymer consisting of hard spheres, the integra- 
tion over Sp (r) yields the following Hamiltonian for the 
coarse-grained field in the grand-canonical ensemble |^ : 

H [{uk} ; {wj] ~ ^ [-/i + A^,,(wj)] rii - e ^ riirij. 



<hj> 



(1) 



The quantity p is the solvent chemical potential, e is 
the nearest-neighbor coupling energy parameter of the 
solvent, the sum over includes all nearest neighbor 
cells, Vi denotes the volume excluded by the solute in cell 
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and 



A^„,(ui) « CVi 



(2) 



with c = 2.67x 10^ J/m'^, is the reversible work to accom- 
modate that excluded volume in the liquid. Equations (|^) 
and (2) are approximations to the full results derived in 
Ref. |8|. Specifically, in comparison with Eqs.(19) or (23) 
of Ref. 1^, we see that Eq.(|l|) neglects relatively small 
inter-cell interactions. Further, Eq.(^) is a reasonable nu- 
merical approximation to the general result for An^^(vi) 
that follows from Gaussian statistics for Sp{r) 0]. The 
specific value of c is chosen so that A/iex(47ri?'^/3) yields 
the excess chemical potential for a hard sphere or bubble 
that excludes solvent from a spherical volume of radius 
R =0.5 nm 

The predictions of the model for excess chemical po- 
tential or solvation free energy of a hard sphere over the 
range of all possible radii is shown in Fig. It is seen 
that the agreement between the results of the model and 
those of an atomistic simulation is reasonable, as the 
model captures the appropriate free energy scaling at 
small and large lengths and gives a reasonable predic- 
tion of the crossover length between the two regimes. 
For closer agreement between the statistical field model 
and the atomistic model in the crossover regime, terms 
neglected in Eq.(|l|) should be included, as demonstrated 
in Ref. g. 

The crossover occurring around i? ~ 1 nm is a collec- 
tive effect. It is induced by the entropic cost of constrain- 
ing small length scale fiuctuations near a hydrophobic 
solute. For a given cell i, the cost is AiJ,^^{vi). When 
the Wi's form a connected volume of large enough extent, 
the net entropic cost is larger than the energetic cost for 
forming an interface that surrounds the volume. At this 
point, and for larger volumes, it becomes energetically 
preferable for m — for all cells i occupied by or adja- 
cent to the large solute, and the free energetics is then 
dominated by the energy of the interface thus formed. 
Such interfaces are stable only in the neighborhood of 
vapor-liquid phase equilibrium. Therefore, crossover be- 
havior can occur only when fi — ficocx ^ ksT, where fc^T 
is Boltzmann's constant times temperature (the thermal 
energy) and /Zcoex is the chemical potential at phase co- 
existence between liquid and vapor. Further, since a 
pressure-volume term generally contributes to the free 
energy of a large void, crossover behavior occurs only 
when pressure is low enough to make this term similarly 
small. We will see below that the crossover in the scal- 
ing behavior of the free energy is closely related to the 
nucleation of oil-water demixing and the mechanism for 
the collapse of a hydrophobic chain. 



DYNAMICS OF THE SOLVENT AND CHAIN 

The excluded volumes Vi are dynamical, changing as 
the chain configuration changes with time. The chain 
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FIG. 1. The solvation energy per unit surface area for 
a liard sphere or bubble as a function of R, where R is 
the radius of the spherical volume from which the solvent, 
water, is excluded. The circles are the results of com- 
puter simulation for the SPC/E atomistic model of water 
[ p^ . The line is the prediction of our model computed by 
Monte Carlo umbrella sampling |Q, using the weight func- 
tional exp J/ [{rife} , {ui}]} , where l3~^ = ksT, and with 
= 47r_R''/3.The parameters of the model have been 
chosen so as to match the isothermal compressibility, the 
chemical potential and the surface tension of water at room 
temperature and 1 atm pressure. This yields I = 0.21nm, 
e = 3.74kJ/mol (= l.SlfcsT at room temperature), and 
M — /icoox = 0.557J/mol (= 2.25 x 10~'*fcsT at room tem- 
perature), where /icoex is the chemical potential of the solvent 
when it is in coexistence with the vapor at room temperature. 
The comparison of — /icoex with e of shows that the bulk 
liquid is extremely close to phase coexistence. The horizontal 
dashed line is placed at the value of the liquid-vapor surface 
tension and coincides with the large R asymptotic value of 
A/icx/47r_R^ for both the model and the atomistic simulation 
(the surface tension of SPC/E water model is very close to 
that of real water ji^] ) . 



need not move slowly on the time scale for which dis- 
turbances in the solvent may relax. As such, reversible 
work surfaces, like the solvation free energy graphed in 
Fig. 1^, are not entirely pertinent, and the dynamics of 
the liquid and chain should be studied together. For this 
purpose, we have constructed a stochastic dynamics that 
satisfies the basic requirements of time reversibility and 
the preservation of an equilibrium distribution. With this 
dynamics, solute particles move through space with con- 
tinuous variation of coordinates while the coarse grained 
density evolves simultaneously through its discrete con- 
figurations. Since 6p (r) has been integrated out (its ef- 
fects appear only implicitly) the trajectories can be true 
to nature only for times greater than those required to 
relax Sp{r). This time is of the order of 1 ps On 
this time scale, the appropriate stochastic dynamics is 
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non-inertial. 

Trajectories of our model are advanced as follows: The 
polymer, composed of = 12 spheres is moved for Mg 
time steps from a configuration {ro,} to a new configu- 
ration {r^} according to propagation rules of Langevin 
dynamics, in the field of constant coarse grained density 
variables, {ni}. At the end of those steps, the polymer 
is held fixed and the coarse grained density is moved to 
a new configuration {n^} by carrying out a full sweep 
through all cells, applying Glauber Monte Carlo dynam- 
ics to each cell with the Hamiltonian H [{n^} , 
The net procedure taking the system from {ra,ni} to 
{r^,nj} is then repeated over and over again. The tra- 
jectories thus formed are reversible, preserve the norm of 
the configurational distribution function, and obey de- 
tailed balance. 

For one of the Alg time steps between solvent sweeps, 
carrying the monomers from their configuration at time 
t to that at time t + 6ts, the Langevin propagation cor- 
responds to 

r„(i + Sts) = ra{t) + Sts-f-'^Fa{{r^{t)} , {n,}), (3) 

where Fa{{r^{t)} , {ui}) is the force acting on monomer 
a, and 7 is the friction constant. The force contains a 
random part, (5F, that is independent of configurational 
variables. It is the dynamical remnant of the small length 
scale field, Sp{r). We take it to be isotropic, to have a 
correlation time of zero, and to obey Gaussian statistics 

with zero mean and variance ^|(5F|^^ = SkBTj. We esti- 
mate the value of 7 by associating the diffusion constant 
of a single bead moving in water, Dg = ksT/^, to its 
Stokes estimate which gives 7 = 8.39 x IO~^^kg/s. The 
remaining contributions to the force is not random, and 
is given by — [y ({r^}) -I- WiJii]. The time step 
is Sts — 1.4 X 10^^** s, which is small enough that (|) 
preserves the canonical distribution. 

A physically meaningful value of Mg coincides with 
Ms ~ 5ti/6ts, where 6ti is the physical time associated 
with a sweep through the coarse grained density of the 
solvent. This solvent time should coincide with the corre- 
lation time for a density fluctuation of length scale I, i.e, 
6ti — 1/[D{2tt/1)'^], where D is the self diffusion constant 
of liquid water. This relationship gives 6ti — 5.0 x 10" ^^s, 
and thus Mg — 36. 

With these dynamical rules, we show below that the 
half-life of the extended chain we consider is approx- 
imately 10"^ s at room temperature, while the tran- 
sient time to move off a threshold between extended coil 
and compact globule states is of the order of 10~^° s. 
Thus, collapse transitions in this system are rare events 
. Harvesting a representative ensemble of these events 
without prior knowledge of transition states or mecha- 
nism is accomplished with transition path sampling [pT| . 
To apply this technique, we use the number of vapor 
cells, Si(l — Ui), as the order parameter to distinguish 
the solvated coil and globule states. 



RESULTS 

Figure ^ illustrates a 1.5 ns trajectory in which a chain 
of 12 hydrophobic units collapses from an extended coil 
to a compact globule. The volume occupied by each sep- 
arate unit is typical of that for amino acids in solution. In 
vacuum, this particular chain will always remain in a coil 
state since the global minimum in its intrachain potential 
energy function corresponds to the fully extended confor- 
mation. Nevertheless, the collapsed chain is very much 
the stable state of the solvated chain, as is evident from 
the free energy surface shown in Fig. ^ While rare, the 
trajectory illustrated in Fig. |^ is a representative exam- 
ple of the sub-ensemble of all 1.5 ns trajectories that do 
exhibit the collapse transition at temperature T w 300 K. 
The polymer is initially in the extended coil state, after 
which a spontaneous fluctuation in chain conformation 
collapses a section of the polymer. The collapsed sec- 
tion forms a sufficiently large hydrophobic cluster that 
the formation of a vapor bubble is induced. Eventually, 
this vapor bubble grows and drives all the units of the 
polymer together. 

Visual inspections of this and similar dynamical trajec- 
tories suggest that the mechanism of the collapse tran- 
sition arises from an interplay between the size of the 
polymer and the formation of a vapor bubble. For this 
reason, we have mapped out the free-energy landscape 
in terms of variables manifesting the polymer size and 
the bubble size. The first of these variables, Rg, is the 
squared radius of gyration of the polymer. The other, 
U, is the volume of the largest vapor bubble, in units of 
l^. In order to identify the bubbles in the system, we 
use the criterion that any two vapor cells (i.e., cells for 
which rii = 0) belong to the same bubble whenever they 
are nearest neighbor cells. 

The contour plot for the room temperature free energy 
landscape. Fig. ^, shows that the path of lowest free en- 
ergy from the coil to the collapsed globule is one where, 
initially, the radius of gyration decreases, while the size of 
the largest bubble is still essentially zero. In this regime, 
the solvent still wets the polymer (i.e., rii is mostly 1 for 
cells occupied by the solute), and the free energy hardly 
changes. When the radius of gyration becomes small 
enough, however, a bubble starts to grow. Here, the 
free energy increases sharply by about dksT where it 
reaches a saddle point at {U,B?g) = (98, 23.5P). Beyond 
that saddle point, the bubble grows spontaneously to a 
size that eventually envelopes the fully collapsed glob- 
ule. The height of the barrier is only weakly dependent 
upon chain stiffness since it is due primarily to solvent 
reorganization. 
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FIG. 2. Four configurations from a trajectory wliere a 12 
unit hydropiiobic drain in water goes from the coil to the glob- 
ule state, (a) shows a configuration from the equilibrated coil. 
The chain remained in configurations like that throughout a 
10 ns run at room temperature (T = 0.663e). On a much 
longer time scale, about 10~^ s, the chain typically does ex- 
hibit a transition from the coil to globule. Such events are 
found with transition path sampling, equilibrating from an 
initial high temperature (T — 0.74e) 10 ns trajectory that 
exhibited the transition spontaneously. Three configurations 
exhibiting that transition, covering a stretch of 1.5 ns are 
shown in (b), (c) and (d), with that in (c) being a configura- 
tion from the transition state surface. The transparent cubes 
denote the vapor cells. Those seen far from the chain are typi- 
cal spontaneous density fiuctuations in bulk liquid water. The 
size of the simulation box is 397 nm"^, corresponding to 42,875 
cells. The parameters characterizing the energy of the solvent 
are given in the caption to Fig.l. The intra-chain potential, 
1/ (ri, r2, ri2), is a function of the positions of the centers of 
each of the 12 hydrophobic spheres (the red particles in the 
figure). It contains three parts: (1) steep (essentially hard 
sphere) repulsions between solute particles such that their 
interparticle separations are larger than a = 0.72 nm; (2) 
stiff harmonic potentials bonding adjacent particles in the hy- 
drophobic chain, |fcs(cr — Ir^+i — fal)^, with ks = 14.1J/m^; 
(3) a bending potential favoring an extended chain, ^kgOa^ , 
where da is the angle between (ra+2 — Va+i) and (r^+i — r^), 
and kg = 1.85 x 10~^°J/rad'^. The volumes Vi excluded from 
water by the chain are dynamic as they change with changing 
chain configuration, i.e., Vi — Vi {{ra}). Specifically, these 
volumes are computed by assuming water molecules have van 
der Waals radii equal to 0.14nm, and that the diameter of 
each hydrophobic unit is cr = 0.72 nm.. That is, points in 
the excluded volume, r, are those in the union of all volumes 
inscribed by |r — Tq | < 0.5 nm, a = 1, 2, A^s. 

A close examination of the free-energy landscape re- 
veals a small barrier of height 2k bT at {U,R^) = 
(6.5, 70P). This feature, not visible on the scale of the 
graph plotted in Fig. ||, separates the coil from a more 
compact metastable intermediate. The presence of these 
two states arises from a competition between the entropy 
of the chain, which favors the coil state, and weak deple- 
tion forces, which favor a more compact state. Weak de- 
pletion forces are caused by the reduction in the volume 
from which the solutes exclude the solvent, when solutes 
come together while still wet. The attraction between 
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FIG. 3. Contour plot of the free energy landscape for 
the collapse of our hydrophobic polymer, computed by 
Monte Carlo umbrella sampling |^|, using the weight func- 
tional exp{—f3H [{rifc} , {vi}] — 13V [{ra}]}, where the intra- 
chain potential, V [{re}], and the excluded volumes, Vi [{r^}], 
are determined as described in the legend to Fig. 2. The 
curves of constant free energy are drawn as a function of the 
squared radius of gyration and the size of the largest bubble 
in the system (see text). Neighboring lines are separated by 
2.5kBT. Superimposed is a scatter plot (in black) of the har- 
vested 150 ps trajectories going from the coil to the globule 
state. The transition states are indicated in green. The har- 
vesting was performed with transition path sampling, making 
8,400 moves in trajectory space, of which 75% were shooting 
and 25% were shifting [0. We find that the plateau regime of 
the flux correlation function is reached after 50-70 ps [ ^"gIJi^ 
implying that the typical commitment time for trajectories 
to pass over the barrier is of the order of 0.1ns. Given this 
time and the fact that the figure shows the free energy bar- 
rier separating the extended coil and compact globule states 
to be about 9fcsT, the half life of the extended chain is about 
0.1 ns X exp(9) ^ 10"^ s. 

two small hydrophobic objects predominantly arises from 
this effect |^,^,^ , though its full description (with the 
characteristic oscillations in potentials of mean force) re- 
quires a more accurate evaluation of A/iox (fi) than Eq. ||. 
Whether evaluated to high accuracy or approximately, 
this driving force is relatively small, and it is not respon- 
sible for spontaneous assembly of hydrophobic units [ pO[ . 
The free energy difference is SO/c^T between the coil and 
the fully collapsed globule, whereas it is only a few ksT 
between the coil and the intermediate state. The strong 
driving force for the collapse of the polymer comes from 
the emptying of cells (i.e., solvent cavitation or drying) 
and thus the demixing of the hydrophobic units and wa- 
ter. The resulting stabilization follows from the fact that 
the solvent is close to liquid-vapor equilibrium so that 
the large length scale interfacial free energy of the cav- 
ity is far lower than the small length scale entropic cost 
of maintaining a wet state, Si Afj,ex{vi)- Thus, nucle- 
ating the collapsed hydrophobic chain involves the same 
physical effect that is responsible for the crossover phe- 
nomenon discussed above in reference to Fig. |l|. 
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To judge the extent to which i?^ and U describe the 
pertinent dynamics of the collapse, and are thus good 
reaction coordinates, we have performed extensive tran- 
sition path sampling of an ensemble of trajectories, each 
of length 150 ps. Points taken from these transition tra- 
jectories are shown in juxtaposition with the free energy 
surface in Fig||. While ten times shorter than the one il- 
lustrated in Fig.|[ the 150 ps trajectories are sufficiently 
long to capture the mechanism of the collapse because 
150 ps is significantly longer than the time for a trajec- 
tory starting at the dynamical bottleneck to commit to a 
basin of attraction. Figure ^ shows that the fiow through 
the bottleneck in this system is primarily due to the sol- 
vent, and the solvent in this transition state regime moves 
on the time scale of picoseconds. To identify the bottle- 
neck (i.e., transition state surface), we have located the 
configurations on each trajectory from which newly initi- 
ated trajectories have equal probability of ending in the 
coil or globule states. These points in configuration space 
are members of the transition state ensemble [^^^ ■ We 
project them onto Fig. ^ where it is seen that points in the 
transition state ensemble are indeed close to the saddle 
point in the free energy surface. Thus, U and Rg are the 
predominant reaction coordinates for this system. The 
transition state ensemble is slightly tilted in the [U, Rg)- 
plane, showing that the larger the polymer, as measured 
by its radius of gyration, the larger the size of the critical 
vapor bubble. However, the scatter of the transition state 
ensemble from a line in that plane is notable, indicating 
that at least one other variable in addition to U and R^ 
plays a quantitative role in the reaction coordinate. 

The transition paths shown in Fig. ^ differ from the 
lowest free energy path, and the differences are generally 
larger than fc^T. This behavior shows that the polymer 
and the solvent move on significantly different time scales 
when passing through and moving beyond the dynami- 
cal bottleneck. When a vapor bubble is nucleated, the 
polymer does not respond on the time scale at which the 
bubble is formed. In fact, by taking an artificially large 
value for Mg (3,600) to accelerate the dynamics of the 
chain relative to that of the solvent by a factor of 100, 
the projected dynamical paths closely follow the lowest 
free energy path. 

DISCUSSION 

Trajectories of the hydrophobic collapse are generally 
parallel to U as they pass over the transition state. This 
finding demonstrates unambiguously that the dynamics 
of collapse of a hydrophobic polymer in water is domi- 
nated by the dynamics of water, specifically the collective 
emptying of regions of space near a nucleating cluster 
of hydrophobic species. Further, both the solvent and 
chain remain out of equilibrium throughout the collapse 
transition. Real hydrophobic molecules have some affin- 
ity to water, and vice versa, on account of ubiquitous 



van der Waals interactions. The ideal hydrophobic chain 
considered here, however, has no such property as it is 
composed simply of hard spheres. Nevertheless, we ex- 
pect our conclusions to remain valid for more realistic 
models of hydrophobic species because the underlying 
mechanism of nucleation is largely unaffected by van der 
Waals interactions Our findings would seem also 

pertinent to the mechanism of biological assembly, such 
as protein folding, but to demonstrate so with simulation 
will require an analogous simulation study of a protein- 
like chain. In the past, a statistically meaningful study 
of such folding with explicit account of solvent dynamics 
seemed impractical. The current analysis, however, sug- 
gests that such studies will become feasible by using a 
statistical field model like that employed herein. 

Appropriate models for hydrophobic collapse must ac- 
count for both the physics of small length scale fluctua- 
tions and the physics of phase equilibria and interfaces. 
This accounting is difficult to accomplish without an ex- 
plicit treatment of the solvent, as we do here. Often, 
solvation is treated with implicit solvent models in which 
all solvent degrees of freedom are averaged out in some 
approximate fashion, assuming all these degrees of free- 
dom relax quickly. Some of these models mimic the ef- 
fects of the solvent through estimates of the pair poten- 
tial of mean force between solute units. This approach 
can be correct for describing small length scale effects, 
and it therefore can be correct in situations where the 
solute remains wet. But for collections of hydrophobic 
solutes large enough to induce drying or cavitation, an 
implicit model would need to describe correlations of es- 
sentially all orders, and not simply pair correlations. For 
this reason, a recent study of hydrophobic effects in pro- 
tein folding requires further analysis . Other implicit 
solvent models are based upon estimates of exposed sur- 
face area |2^,Q. This approach can be valid for large 
enough solutes provided the time scales of interest are 
longer than nucleation times (10~^ s in the system stud- 
ied here). But applying it to small unassembled solutes 
is incorrect because the free energy in this regime scales 
with excluded volume and not surface area. 

The connection we stress between hydrophobic col- 
lapse and phase equilibria provides a simple explana- 
tion for both cold denaturation and pressure denatura- 
tion of proteins. Namely, the tendency for drying and its 
concomitant stabilization of hydrophobic clustering de- 
creases with decreasing proximity of liquid- vapor equilib- 
rium. Thus, the lowering temperature and the increasing 
of pressure destabilize hydrophobic collapse because both 
these actions move ambient water away from its phase 
equilibrium with vapor. The proximity of phase equi- 
libria is also of relevance to free energy barriers for hy- 
drophobic collapse, like that illustrated in Fig. ^. These 
barriers will increase, for example, with increasing exter- 
nal pressure. Thus, the most direct measure of the im- 
portance of hydrophobic collapse in protein folding may 
therefore come from studying the kinetic effects of chang- 
ing pressure. 
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